Multiple Vortex Cores in 2D Electronic Systems with Proximity Induced 

Superconductivity 
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The structure of a proximity induced vortex core in a two-dimensional (2D) metallic layer covering 
a superconducting half-space is calculated. We predict formation of a multiple vortex core charac- 
terized by two-scale behavior of the local density of states (LDOS). For coherent tunneling between 
the 2D layer and the bulk superconductor, the spectrum has two subgap branches while for incoher- 
ent tunneling only one of them remains. The resulting splitting of the zero-bias anomaly and the 
multiple peak structure in the LDOS should be visible in the tunneling spectroscopy experiments. 

PACS numbers: 73.22.-f; 74.45. +c; 74.78.-w 



Experimental study of subgap quasiparticle states in 
a superconductor (SC) placed in a magnetic field pro- 
vides a unique tool for probing the internal structure of 
Cooper pairs. These states bound to the vortex cores are 
strongly affected by the superconducting gap anisotropy 
in the quasimomentum space (see [lj and refs. therein) 
and by the number of the order parameter components^. 
Direct information of the spectrum and of the wave func- 
tions of such excitations can be obtained by scanning 
tunneling microscopy/spectroscopy (STM/STS) which 
probes the energy and spatial dependencies of the local 
density of states (LDOS)^. However, the existing exper- 
imental data often provide rather controversial informa- 
tion on the so-called zero bias anomaly known to be a 
fingerprint of the Caroli-dc Gennes Matricon (CdGM) 
states within the vortex core^. An obvious reason for 
such ambiguity can be a defect surface layer which masks 
the bulk quasiparticle states. For example, an energy gap 
in such (possibly nonsuperconducting) layer can appear 
due to proximity to the superconductor (see, e.g., [f|). 
The vortex states in the systems with proximity induced 
superconductivity have been recently studied using vari- 
ous phenomenological approaches aimed to describe the 
hybrid structures consisting of graphene layers coupled 
to superconducting electrodes^—. 

Here we propose the microscopic description of the vor- 
tex core states in two-dimensional (2D) electronic sys- 
tems with proximity induced superconductivity and an- 
alyze the masking effects of a thin surface layer on the 
STM/STS data. Based on the general approach devel- 
oped in Ref.— for proximity induced superconductivity 
we formulate two models which describe the electron 
transfer between a 2D system and a bulk superconduc- 
tor: (i) coherent momentum-conserving tunneling model 
and (ii) incoherent tunneling model that accounts for dis- 
order and corresponding breakdown of momentum con- 
servation for tunneling quasiparticles. Within both these 
models, proximity to a superconductor induces supercon- 
ducting correlations in the 2D layer and leads to forma- 
tion of an energy gap with a magnitude depending 



on the tunneling rate T— ~— : ~ T for T <C A, where 
A is the gap in the superconducting electrode. The hall- 
mark of the induced gap is that it does not have a sep- 
arate critical temperature but rather vanishes together 
with the bulk gap. Quite naturally, the spatial behavior 
of quasiparticle wave functions in the 2D layer is deter- 
mined by two length scales: (i) the coherence length, 
£s = HVp/A for clean or £5 = yf hD$ / A for dirty limit, 
characterizing the bulk electrode and (ii) the 2D coher- 
ence length, &d = KV 2D /A 2D or £_ 2D = ^/hD 2D /A 2D , 
where Vp, V%r> and D$, D 2 d are the Fermi velocities 
and diffusion constants in the bulk and in the 2D layer, 
respectively. Since A 2 r> < A the coherence length £ 2 d 
usually is longer than £5. 

We show that the proximity induced vortex in a ballis- 
tic 2D layer has a "multiple core" structure characterized 
by the two length scales, £g and £2,0- Such a two-scale 
feature did not appear in the preceding theoretical works 
where proximity vortex states have been induced by a pri- 
mary vortex pinned at a large-size hollow cylinde r 14 ! 15 . 
We calculate the energy spectrum of core excitations for 
both coherent and incoherent tunneling. For coherent 
tunneling, the spectrum of quasiparticles bound to the 
multiple core consists of two anomalous branches as func- 
tions of the impact parameter b. One branch, £\{b), qual- 
itatively follows the usual CdGM anomalous spectrum 
eo(b) of the primary vortex; it extends above the induced 
gap where it turns into a scattering resonance. The other 
branch, e 2 (b), lies below the induced gap and resembles 
the CdGM anomalous spectrum for a vortex with a much 
larger core radius ~ £2,0 ■ This branch has a much slower 
dependence on the impact parameter and reaches the in- 
duced gap for trajectories that completely miss the core 
of the primary vortex. 

We demonstrate that the structure of the multiple core 
is strongly affected by disorder, i.e., by impurity scatter- 
ing inside the bulk electrode or inside the 2D layer, as well 
as by the barrier disorder. In our incoherent tunneling 
model, the latter is accounted for by ensemble averag- 
ing over various realizations of disorder in the barrier, 
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which results in the suppression of influence of the pri- 
mary CdGM spectral branch on the spectral characteris- 
tics of the 2D layer. The lower anomalous branch 62(b) 
survives the destructive influence of the barrier disorder, 
though being shifted and broadened due to the momen- 
tum uncertainty. The impurity scattering inside the bulk 
superconductor and inside the 2D layer causes further 
smearing of the spectral characteristics of the core states 
which approach the usual LDOS for dirty superconduc- 
tors scaled with the corresponding coherence lengths £213 ■ 
As a result, the spatial and energy dependence of the 
LDOS inside the multiple core reveals a rich behavior 
dependent on the above spectral properties. The LDOS 
contribution from the larger core region ~ ^2D behaves 
similarly to the standard vortex LDOS with the corre- 
sponding gap and coherence length. Effects of primary 
vortex core on the spatial LDOS pattern in the 2D layer 
can be seen as a narrow peak which strongly depends on 
the degree of disorder. 

Model. Hereafter we consider a superconducting half 
space coupled by quasiparticle tunneling to a 2D cover- 
ing normal metal layer. We start from the quasiclassical 
Eilenberger equations for the retarded or advanced Green 
function in the 2D layer which are easily derived from the 
results of Ref.— (see Appendix |X] for details) : 

- ihv 2D V5(p 2 D, r) - e [f 3 g(p 2D , r) - g{p2D, r)f 3 ] 

- [St<?(P2D, r) - g{p 2 D, r)S T ] = , (1) 

where P2D and v 2 r> = <9e2,o(p)/<9p are the 2D layer 
Fermi momentum and velocity. The Pauli matrices n, 
T2, T3, the Green function 
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and the self energy are matrices in the Nambu space. 

Coherent tunneling conserves in-plane momentum 
P2D- Therefore, the self energy St has the form 

S T (P2D,r) = % — [3s(p+;r,0)+$ s (p_;r,0)], (2) 
where 3D momentum p± = (p2T>, ±Pz) lies on the Fermi 



surface of the bulk SC, p\ D + p z 



p\- 



Coherent tun- 



neling is impossible if the Fermi momentum in the 2D 
layer is larger than that in 3D, p 2 D > Pf- For incoherent 
tunneling, the in-plane momentum is not conserved. All 
momentum directions participate in tunneling, thus 



t T (r)=iT(gs(p F ;r,0)) 



(3) 



The angular brackets denote the averaging over direc- 
tions of the 3D Fermi momentum p^ . The tunneling rate 
r ~ t 2 1 E2F can be expressed-^ in terms of the normal- 
state tunnel conductance G — 1/RS per unit contact 
area, L = G / (AttGo^) ~ E 2 fRo/R 7 with the conduc- 
tance quantum Gq — e 2 /irh and the normal 2D density 
of states (DOS) v 2 = m/2-Kh 2 . Therefore T/E 2F < 1 
if the total tunnel resistance R is much larger than the 



Sharvin resistance Ro = (NGq) -1 for an ideal ./V-mode 
contact with the contact area Si£. Nevertheless, there is 
a room for the condition T ~ A to be fulfilled even for a 
large contact resistance R ^> Rq. 

Here we restrict ourselves to the limit of low tunnel- 
ing rate T <C A which leads to a small induced gapii 
A 2 t> = F and long coherence length £ 2 r> 3> £s- We con- 
sider an isolated vortex line oriented along the z axis 
perpendicular to the SC/2D interface and choose the 
gap function inside the bulk superconductor in the form 
A = Ao(/o)e 1 *, where (p,cj>,z) are the cylindrical coor- 
dinates; Aq(p) approaches the bulk value Aoo far from 
the vortex core. The self energies in the 2D layer are 
given by Eqs. © and ([3]). They have parts with sharp 
peaks localized at small distances p ~ £s and the adia- 
batic "vortex potential" part IS.2D ~ Te 1 ^ which defines 
the large scale behavior of the 2D layer Green functions. 

Multiple core. Clean limit with coherent tunneling. To 
elucidate the basic features of the multi-scale vortex core 
in the 2D layer we consider first an idealized picture with- 
out any disorder assuming specular electron reflection at 
the surface of the bulk SC. 

For the low-energy limit e <C A^ one can find the 
induced vortex potential at large distances p> £,g : 
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(4) 



The quasiparticles propagating along the trajectories 
that miss the primary vortex core (b > £s) are affected 
only by this long-distance (£2.0 3> £s) P ar t of the in- 
duced gap potential and the corresponding solutions for 
the Green functions coincide with the standard CdGM 
ones for the gap value replaced with T. A quasiclassical 
trajectory can be parameterized by its angle a with the 
x axis, the impact parameter b — psin(0 — a) and the 
coordinate s = pcos((f>—a) along the trajectory. The cor- 
responding anomalous spectrum for 2D excitations i a 17 i 18 
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where A = £2d/H- This modified CdGM branch should 
dominate in the local DOS at large distances p 3> <ts- 

Trajectory with a small impact parameter b < £s can 
be divided into the long-distance part going far from the 
primary vortex core, and the region inside the core. Far 
from the core, the solution is found using the vortex po- 
tentials Eq. In the core region one should take into 
account the self-energy parts localized within the pri- 
mary vorte x 17 ' 18 . We put St = *L (Fo • r), where f — 
(fie-^ a ,T 2 e-^ Q ,f 3 ), the vector F = (-Cs,9 s ,g S ) is 
normalized by Fg = 1 and has the components 
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2Q [e - eo ± iS\ ' 

e (b)=bQ- 1 I [Ac/ p]e~ K ^ds , (7) 
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Q = / e- K ^ ds ; K(s) = — / A (p') dp' . (8) 
Jo n^H J\b\ 



3 




FIG. 1: (Color online) Two-scale behavior of the spectrum, 
Eq. ([9]), for coherent tunneling. The spectrum has two local- 
ized branches, £i(6) and £2(6), for e < T. The branch £2(6) 
has a scale £20, it saturates at e = V for b ^> &d- £i(&) has 
a scale £s- For 5fte > T it transforms into scattering reso- 
nances, b* is defined as £i(6*) = T — 0, while 6' corresponds 
to SRei(fe') = T + 0. 
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FIG. 2: (Color online) The local DOS in logarithmic scale for 
coherent (right panel) and incoherent (left panel) tunneling 
in the clean limit for different distances p from the vortex 
center. The peaks in LDOS exist up to distances ~ £20- 
Here A/T = 5, q = 1. 



Here V|| is the 3D Fermi velocity projection onto the (x, y) 
plane. The upper (lower) sign of an infinitely small 6 > 
refers to the retarded (advanced) function. Matching the 
2D Green functions through the primary core region of 
rapidly changing self-energy potentials (see Appendix |C|) 
we find both the spectrum 



q ^[e-e 2 (b)][e-e Q (b)} 



£ 2 (&)] 2 =0(9) 



and the Green functions for trajectories with 6 <C £213 • 
Here q = v\\/V2D- For b < £s, the cut-off parameter in 
Eq. (0 should be replaced with A = &D/S.S- 

The resulting two-scale behavior of the spectral 
branches is illustrated in Fig. Q] The complex- valued en- 
ergy branches satisfy the symmetry condition: £1,2(^6) = 
— e 12 (&). There are two real-valued energy branches in 
the range |e| < T crossing zero of energy as functions of 



the impact parameter and one complex- valued branch in 
the range T < |e| < Aoo. The lowest-energy branch e 2 (b) 
as a function of the impact parameter has a character- 
istic scale ^2D- For b < it is determined by Eq. ([5]) 
with the proper cut-off parameter A as discussed above. 
On the other hand, it saturates at e = T for b >> 
The branch £1(6) has a scale £s: At low energies it goes 
slightly below the CdGM spectrum eo(b) in the bulk SC, 
e i(b) = (1 + q/2)~ 1 eo(b). Above the induced gap T the 
spectrum transforms into a scattering resonance due to 
the decay into the delocalized modes propagating in the 
2D layer: £1(6) = e (6) - iT for |e| > T. Since Eq. © 
determines a pole of the retarded Green function in the 
lower half-plane of complex £, the square root in Eq. © 
should be analytically continued under the cut extending 
from —00 to — r and from T to +00. As a result, £1(6) 
has a discontinuity at e± = T. 

Two branches appear because the system under con- 
sideration consists of two sub-systems^, the bulk SC 
and the 2D proximity layer, each with its own anoma- 
lous branch. The branch £1(6) is the proximity image of 
the bulk spectrum £q(&) with the spectral weight propor- 
tional to the tunneling probability T. The branch £2(6) 
belongs to the 2D layer itself. We note that the presence 
of two anomalous branches does not contradict to the 
index theorem^. Indeed, its application requires that 
both zero of the quasiclassical Hamiltonian at the Fermi 
surface and its singularity at e — eo(b) are taken into 
account when calculating the topological invariant. As 
a result, the number of anomalous branches is increased 
up to 2 for a single-quantum vortex. 

The multiple-branch spectrum results in multiple 
peaks in the LDOS energy dependence (right panel in 
Fig. [2]). Here the local DOS is defined by the angle- 
resolved one (normalized by its normal state value) 
N e (s,b) = [g R (s, b) — g A (s, b)]/2 averaged over the tra- 
jectory direction. The multiple peak structure appears 
to be most pronounced deeply inside the primary core 
region (at distances p < Q/^d when t\ < T) illustrat- 
ing, thus, the two-scale structure of the vortex core. 

The number of LDOS peaks at a certain distance p 
from the vortex center is determined by the number of 
spectral branches at b ~ p. The spectrum discontinu- 
ity at the induced gap T causes the appearance of three 
LDOS peaks in the range of distances, corresponding to 
b' <b<b* (see the plot for p = 0.3£ s in Fig.|2|). The nu- 
merical LDOS patterns have been obtained by the sub- 
sequent solving of two sets of Eilenberger equations in 
Riccati parametrization 21 : first, we calculated the Green 
functions for the bulk superconductor with the model or- 
der parameter profile A (p) = A^p/ vf+^s an d, sec- 
ond, we have found the solution of Eq. (JXJ) for a 2D layer 
with the induced potentials defined by the Eq. ©. 

Multiple core. Clean limit with incoherent tunneling. 
We proceed our study with the consideration of disorder 
effects and introduce first the momentum scattering dur- 
ing the tunneling process described within the incoherent 
tunneling model. Considering the tunneling as a pertur- 
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bation one can assume a specular quasiparticle scattering 
at the interface and, thus, use the results of the previous 
section for the Green functions. The self-energy poten- 
tials in this case can be obtained by averaging of Eqs. 
([6][8]) over the trajectory direction: Ej- = *r(.9s)- This 
averaging does not affect, of course, the induced gap func- 
tion Q outside the primary vortex core and, thus, the 
spectrum e 2 survives the influence of the tunnel barrier 
disorder at least for b > £5. On the contrary, the subgap 
branches localized within the primary vortex core region 
appear to be completely destroyed. Such dramatic con- 
sequence of the momentum scattering is caused by the 
averaging of electronic wave functions with different im- 
pact parameters and consequent loss of any information 
about the CdGM states of the primary vortex. A natural 
consequence of the momentum scattering is the appear- 
ance of a finite broadening of energy levels for trajectories 
with small impact parameters b < £g. 

Following again the matching procedure described 
above (see Appendix [Cjl we find the angle-resolved DOS 
for b < £ s and |e| < T, 



N e (s,b) 

m 

7(6) 



r 7 (6) £ 



[e-e 2 (b)-(3(b)] 2 +7 2 (b) ' 



r 2 (^sign( £ + & ) 



In 



\QQ \Qb + e\ 



(10) 

(11) 
(12) 



The angular brackets denote averaging over the momen- 
tum p z along the vortex axis in bulk SC, il = deo/db. 
The DOS has a peak of the height r/7 at an energy 
e — 62(b) + (3(b) shifted from a standard bound state level 
(see Appendix IC 21 for details). This shift results in split- 
ting of the zero-bias anomaly in LDOS, as is seen from 
our numerical analysis (see the left panel in Fig. [5]). For 
LDOS calculations we use the numerical procedure sim- 
ilar to that used for the coherent limit above with the 
induced potentials averaged over the Fermi surface (as- 
sumed cylindrical) in the bulk. 

Multiple core. Dirty superconductor with clean 2D 
layer. Smearing of the energy dependence of the in- 
duced potentials caused by disorder becomes even 
stronger if the bulk SC has short mean free path: I <C £s- 
In dirty limit, the momentum averaged retarded (ad- 
vanced) Green functions are parameterized as follows: 



9 S \P) 



r 3 sm 



qR(A) 



T2 COS 



(13) 



We put Q R ( A ) = ±Oi + 162- The boundary conditions 
(|4]) require 61 -> tt/2, 2 ->• for p -> 0. At large 
distances sinOi —> 0, tanh02 — e/T for e < Aoc while 
cosOi — > 0, tanh02 — > — Aoo/e for e > A^. Therefore, 
O2 = for e -c Aqo, and the Usadel equation becomes^, 



Dc 



V 2 Oi 



sin(26i) 
2p 2 



- 2A sin6i = 



(14) 




FIG. 3: (Color online) The local DOS in logarithmic scale for 
the dirty limit with the parameters A/T = 5, V2d/Vf = 1 for 
different distances p from the vortex center. 



The solution of Eq. (fLfj) has been found in Refi^ 2 -: the 
function Oi(p) monotonously decays from 7r/2 at the 
origin down to the zero value at p ^> £g. The Green 
functions (|13l) determine the induced vortex potentials 
Et = iFgs- For |e| -C T and b < £s the peak in the en- 
ergy dependence of the angle-resolved DOS is described 
by the Eq.^ with (3 = and 



2r 2 



2D Jq 



sin 61 ds 



(15) 



(see Appendix IC 31 for details). The numerical results 
clearly confirm the existence of one broadened peak in the 
LDOS dependence vs energy: this peak shifts with the 
increasing distance from the vortex center and becomes 
sharper (see Fig. [3]). Our numerical procedure of the 
LDOS calculation in this limit is based on the using of 
a standard relaxation method 23 for solving the Usadel 
equation 24 in the bulk SC and Riccati parametrization 
for Eilenberger equations in the 2D layer. 

Vortex core expansion. Dirty superconductor and 2D 
layer. To complete our analysis we discuss the case of 
strong disorder both in the bulk superconductor and in 
the 2D layer. This limit has been previously studied in 
ReL2£. As before, one can parameterize the Green func- 
tions averaged over the 2D momentum directions in the 
form of Eq. (fTS"]) , where we use ^ for the 2D-layer func- 
tions instead of O. The boundary conditions coincide 
with those for Eq. (IT3"|) where A^ is replaced with T. 
With the self energies from the previous subsection, the 
Usadel equation for the retarded function for e <C Aoo is 



D 



2D 



V 2 f 



sin(2^) 



2p 2 



2rsin(* - 6) - 2iecos* = 0. 

(16) 

O is essentially nonzero only inside the primary core re- 
gion p < £s- The condition £s < &d = \/UD2d]T 
ensures that such short-distance inhomogeneity in the 
induced vortex potentials does not disturb the adiabatic 
solution based on Eq. §4§ (see Appendix [CM]). Thus, 
putting O = in Eq. (I16p we reduce our problem to 
that describing a standard vortex in a dirty supercon- 
ductor with the gap value T. Thus, the full disordered 
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system should reveal the same LDOS patterns as in the 
bulk case, though scaled with the much larger coherence 
length £ 2 d instead of £g. This conclusion is, of course, in 
agreement with numerical calculations in the Refj^ 5 .. 

Conclusion To summarize, we calculate the electronic 
structure of a proximity induced vortex core in a 2D 
metallic layer covering a superconducting half-space. We 
predict formation of a multiple vortex core resulting in a 
two-scale behavior of the LDOS. For coherent tunneling 
between the 2D layer and the bulk superconductor, the 
spectrum has two subgap branches while for incoherent 
tunneling only one of them remains. The splitting of the 
zero-bias anomaly and the multiple peak structure in the 
LDOS should be visible in the tunneling spectroscopy 
experiments. Disorder further smears the multiple peak 
structure inside the double-scale vortex core. When both 
the bulk SC and the 2D layer are in dirty limits, the 2D 
LDOS qualitatively repeats that in the bulk SC scaled 
with the larger coherence length £20- Such expansion 
of the vortex core probably relates to the anomalously 
large vortex images observed in MgB^ and high- T c 
cuprates2£. 
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Appendix A: Eilenberger equations for coherent and 
incoherent models 

We start with the equation for the retarded (advanced) 
Green functions derived in Refj^i 

G- 1 (r 1 )6 2D (r 1 , r 2 , e) -J E T (n, r')G 2D (r', r 2 , e) dV 

= Lrttn - r 2 ) . (Al) 
Here d is the layer thickness, f±, T2, and T3 as well as 

G-iD = ( ~ ^ , G _1 (ri) = e 2D (p) -H- ef 3 , 

are matrices in the Nambu space, £20 (p) is the spectrum 
of the 2D electron system, p = —iKV and r are the 2D 
momentum and coordinate, correspondingly. The self- 
energy takes the form, 

E T (ri,r 2 ) = dt(ri)(5s(ri,*i = 0;r 2 ,z 2 - 0)t(r a ) (A2) 

where t(r) is the tunneling amplitude which is assumed 
real and the Green function Gs of the bulk SC is taken at 
the SC/2D interface z = 0. The above equations can be 
strongly simplified using a standard quasiclassical proce- 
dure which allows us to derive the Eilenberger equations 



for quasiclassical Green function 

<?(P2D,r) = (m)- 1 I d£ 2 G 2D (p,r) 



Here we derive expressions for the self energies flU |3]) 
and the Eilenberger equations ([1]) for different tunneling 
models. 



1. Coherent tunneling 

Let us assume that the in-plane momentum projection 
is conserved during the tunneling process. This amounts 
for a tunneling amplitude t(r) independent of the co- 
ordinate along the SC/2D interface. In 2D momentum 
representation the self energy in Eq. (jAll) becomes 



t T (Pi,p') = dt 2 J G s (pi,p 2 ;p',^)^^ 



We now apply the operators to the Green function from 
the right and subtract this equation from Eq. (|A1[) . Inte- 
grating the result over the energy variable near the Fermi 
surface £2 = £2d(p) — M and using 



HI 



d 2 p' 



— / ^T(Pi,p')G 2D {p\p 2 )——^ = dt 2 g(p 2D ,r) 



(2tt) S 



x / -^-9s{P2D 1 Pz;r,0)TTiS A [es(P2D,Pz) - l4 

we obtain the quasiclassical Eilenberger equation (fTj). 
The quasiclassical Green function gs(p;r, 0) of the bulk 
SC is taken at the SC/2D interface z = 0. We use the 
mixed momentum p-coordinate r representation describ- 
ing the relative and center-of-mass motion of electrons in 
the Cooper pair and put 

Gs{p,Pz;r,z) = gs(p,Pz;r,z)TriS A (£ 3 ) , 
02£>(p,r) = 9(p,r)m5 A (£ 2 ) . 

Here £3 = es(p,Pz) — /x is the normal quasiparticle spec- 
trum in the 3D half-space, gs and g are standard quasi- 
classical Green functions, and <5a(£2,3) is a delta function 
broadened at the gap energy scale A. 

Assuming isotropic Fermi surfaces in both the super- 
conductor and the 2D layer we get the self energy in the 
form of Eq. $2$ with the tunneling rate 

00 

L = dt 2 J 5 A [e s (p2D,Pz) - (4 dp z . 


Provided the 2D Fermi surface is smaller than the 
extremal cross section of the 3D Fermi surface, i.e., 
P2D < Pf the expression for the tunneling rate reads: 
r = dmt 2 /p3 Z . For large 2D Fermi surfaces p 2 o > Pf 
the self energy term vanishes, and the coherent tunnel- 
ing is impossible. The case of close momenta p 2 u — Pf 
deserves special consideration which should take account 
of a finite delta function width: V ~ (it 2 (m/A) 1 / 2 . 
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2. Incoherent tunneling 

We now assume that the tunneling occurs through ran- 
dom centers such that the ensemble average of ampli- 
tudes in Eq. flXjJ) is 

(i(ri)t(r 2 )) = fsaSin - r 2 ) , (A3) 

where s a is the correlated area of the order of atomic 
scale. After averaging the self energy becomes: 

S r (ri,r 2 ) = t 2 ds a G s (r 1 ,r 1 ;0)6(r 1 - r 2 ) 
= t 2 ds a ini< 3 (0) (g s {p; r, 0)) <5(ri - r 2 ) . 

Here 1/3(0) is the normal density of states in the bulk 
material. Angular brackets denote averaging over three- 
dimensional momentum directions. Within the quasi- 
classical approach the resulting self energy is given by 
Eq. ((3|) with the tunneling rate T = iri> 3 (0)ds a t 2 . This 
approximation coincides with that used in RefJ^. 
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FIG. 4: (Color online) The coordinate frame near the multiple 
vortex core. Primary (induced) core is shown by the white 
(gray) circle. The quasiparticle trajectory with an impact 
parameter b (line AB) passes through the point (p, <f>) shown 
by the black dot. 



Appendix B: Induced vortex potentials 

In both tunneling models the Green functions of the 
2D layer satisfy the Eilenberger equations (TTJ): 

-iftv 2B V/-2[e + S 1 ]/ + 2E 2 5 = 0, (Bl) 
ihv 2D Vp - 2 [e + Ei]/t + 2£ 2 . 9 = 0, (B2) 
-ihv 2D Vg + E 2 p -tfj = 0. (B3) 

and the normalization condition g 2 — ff' = 1 with the 
self energy ([U |3J) as effective potentials 




Quasiparticles are conveniently described by the coor- 
dinates along their trajectories (see Fig. [4j . A quasiclas- 
sical trajectory is parameterized by its angle a with the 
x axis, the impact parameter b = psin(</> — a) and the 
coordinate s — pcos(0 — a) along the trajectory. We 
introduce the symmetric and antisymmetric parts of the 
Green function a 17 ' 18 : 

/ = - [C(s) + i0(s)] exp(ia) (B4a) 
f1 = [C(a)-i6(8)]exp(-ia), (B4b) 

where ((s) = C(~ s ): and 0(s) = —6(—s), The normal- 
ization condition requires g 2 + 9 2 + ( 2 = 1. Eilenberger 
equations (|BHIB3|) can be rewritten as follows: 

fiV r 2D^r + 2(e + £i)0-2i 5 £ fl =0, (B5) 
as 

dO 

W 2D — -2(e + ^ 1 )^-2ig^ I -0, (B6) 
as 

hV 2 D?r + 2iCER + 2i6Ei = 0, (B7) 



where 

2E fl = (s 2 e- M + S 2 e IQ ) , (B8a) 
2iE 7 = (s 2 e" iQ - E 2 e iQ ) . (B8b) 

In order to evaluate the induced vortex potentials Ei, 
E 2 and E 2 we start from two important assumptions: 
(i) low interface barrier transparency and (ii) negligible 
effect of the diffusive interface reflection. These assump- 
tions allow us to neglect the effect of tunneling on the 
bulk superconductor characteristics and use the bulk val- 
ues of the quasiclassical Green functions. Restricting our 
consideration to the small energy values e -C Aoo we find 
the large-scale (p 3> £g) self energy Q to be indepen- 
dent of the particular tunneling model and disorder rate 
in the bulk SC: Ei = iTg s « 0, E 2 = iTf s « Te^, i.e. 
E_r « Ts/p, E/ « Tb/p. 

Contrary to the large distance limit the induced vortex 
potentials close to the primary vortex core reveal a very 
peculiar behavior depending on impurity concentration 
and momentum conservation during the tunneling pro- 
cess. In a clean limit of the bulk SC we use the Green 
function parametrization similar to (|B4[) and rewrite the 
Eilenberger equations in the following form: 

Hvn + 2e6 s - 2iA g s s/p = 0, (B9) 

OS 

tio n - 2e( s - 2iA g s b/p = 0. (B10) 

For energies e <C A^, the functions gs and fs, f s are 
large near the vortex. Assuming that £| 3> 9 2 s — 1, we 

have gg^ = iCs^ • The plus sign here is chosen to 
satisfy the condition of vanishing g$ at large distances 
according to Eq. (@J. The solution of Eqs. (jB9l IBlOj) for 
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retarded and advanced Green function s 17 ' 18 
hv\\e~ K 



MA) 



e 



R(A) 



2Q [e±iS- eo] 
2 

hv 



(Bll) 



[ a (e-A b/p){* A) ds' , (B12) 
Jo 



coincides with ([6][8]). These expressions hold as long as 
\(s\ exceeds \b\/p. For s ^> £,s the function (5 assumes 
its asymptotic expression Cs — ~~ b/p corresponding to 
the boundary conditions ((j). 

In the clean limit the vortex potentials induced in the 
2D layer are crucially dependent on tunneling model. As- 
suming specular quasiparticle reflection at the supercon- 
ductor surface we put gs(+p z ) = 9s(~Pz) m Eq. © so 
that the self energy coincides with the Green function in 
the bulk for coherent tunneling Ei = iTgs, E2 = i^fs 
and with its values averaged over the ensemble for the 
incoherent one: £1 = iT (gs), E2 = ^ (fs)- The ensem- 
ble averaging in terms of quasiclassical Green functions is 
equivalent to the averaging over the 3D momentum direc- 
tion. One can separate two terms in the Green function 
expressions: 
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R{A) 



lC R(A) 



-K 



P 



2Q [e - e ] 



± 



-K 



2Q 



S(e-e ) 



(B13) 

The first term has to be taken as a principal value inte- 
gral when calculating the angular averages. The second 
term is proportional to the delta function of energy and 
determines the density of states (DOS) of the vortex core 
states in the bulk SC. Similarly, the anomalous functions 

f R(A) = e i4^R(A) _ 6 R(A)^ b + is y fi ^ (BW) 
pR(A) = e -i^^R(A) + 9 R(A)^ b _ is y p (B15) 

can be separated into the principal value part and the 
delta-functional contribution. 

Performing averaging over the polar 6 p and azimuthal 
a angles we take into account the symmetry of the func- 
tions under the s-inversion transformation. Thus, we find 
the following expressions for the self energy terms: 



Si = 

E 2 e - ^ = £te # = 



-r (Cs(a)) (B16) 
T(9 s (s)s-<:s(s)b)/p. (B17) 



It is convenient to split the off-diagonal induced poten- 
tial into the localized (E 2 oc ) and the long-range parts: 



E 2 e~ 



= r$ 



$(p) = p(7( a )sign( S )/2Q[e-e ]> 



(B18) 
(B19) 



Here we put I(s) = 2 / Q s (e- A b/p)e~ K ^ ds' . The long- 
range function $ can be regarded as an adiabatic induced 
superconducting gap. Hereafter we focus on the evalu- 
ation of the localized part which is most important in 
the primary core region. Averaging over the azimuthal 



trajectory angle a we find: 



3?e' oc = r. 



hvw 



-K 



2Qilp 



x/e 2 - n 2 P 2 



QE 2 oc = 
5REi = 
9Ei 



±r 



/ ehv\\e K 
\2QpriyJri 2 p 2 - e 2 

/ Hv n e- K 



-sign(e)r 



2 „2\ 



\2QyJe 2 - fl-p 



2 „2 



X (e 2 - n 2 P 



= ±r 



-K 



\2Q^ff 2 p r 



=x(^V-e 2 ) 



Here the upper (lower) sign corresponds to a retarded 
(advanced) self energy term, = deo/db, 



x( x ) = 



1, x > 1 
0, x < 1 



is the Heaviside theta-function, and we use the notation 



sin 8pd8p (. . .) 



for the average over the polar angle 9 P of the 3D Fermi 
momentum. Note that our calculations are essentially 
based on the first-order approximation in the small pa- 
rameter b/p. According to Eq. (|B8|) the symmetrical 
E/(— s) = Ej(s) and antisymmetrical £#(— s) = — £#(s) 
parts of the off-diagonal self energy term X^e -1 ^ can 
be rewritten as follows: E^ 
Eze-^b/p. 



E 2 e-^s/p and E/ 



Appendix C: Scale separation inside the multiple 
vortex core. 

In this Appendix we present the details of the analyt- 
ical procedure used to match the solutions of quasiclas- 
sical equations through the primary core region. In a 
clean 2D layer we start our consideration of the induced 
vortex states from the Eilenberger equations (|B5IIB7[) for 
retarded and advanced Green functions. In the low en- 
ergy limit e <C Aoo appropriate boundary conditions far 
from the induced vortex core (p 3> £20) take the form: 



Vr 2 - e 2 



C = 



Vr 2 - e 2 



,9 = 



Vr 2 - e 2 



(CI) 



The self energy terms reveal a quite different behavior in 
the small (p < £s) and large (p ^> £g) distance regions. 
To match the solutions in these domains we introduce 
a certain distance po such that £s <C po <C C2D and 
consider the Green functions in two overlapping spatial 
intervals: (i) p < p -C £ 2 d and (ii) p > £g. 

Outside the primary core region p^- £,s Eqs. (|B5IIB7|) 
for both tunneling models and arbitrary disorder rate 
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0. 


(C2) 


0, 


(C3) 


. 


(C4) 



inside the superconductor take the form: 

hV 2D ^- + 2e9 - 2igTs/p = 
as 

W 2D ^--2eQ-2igTb/p = 
as 

W 2 J§- + 2i9Tb/p + 2iCTs/p = 
as 

The functions g and £ are even in s while 9 is odd, 
so we can consider only positive s values. We ob- 
tain the solution of the above equations using the first 
order perturbation theory in the impact parameter b: 
w(s) = wq(s) + wi(s), where w(s) = (£,9,g) T . This 
approximation holds for |6| <C £ 2 d- The zero order in b 
solution reads 



w (s) 
where 



1 



u±(s) 



F«o(s) 



c 



Vr 2 - e 2 



u-(s) , (C5) 





Vr 2 - e 2 


■( 


±e 















— it 





„±As 



nv 2D 

This solution satisfies the boundary conditions g — 
-ie/VT 2 - e 2 , C = and 6 = T/^JT 2 - e 2 for s oo 
and e 2 < r 2 . The first order correction wi can be writ- 
ten as 

- )M Cq(s) , C+(s) . C_(s) . 



2 _ ,2 



where 



e- Vr 2 - e 2 Vr 2 - e 2 ' 



(C6) 



6dC (s) = 2C6/ e" As — , (C7) 

P 



t 2D c + ( S )=-b r e - As -, 



6dCL(s) = -& / e 



P 



(C8) 
(C9) 



The lower limit of integration s c in C_ has to be taken 
s c ~ £s for trajectories that go through the primary vor- 
tex core, b < £s, so that the logarithmic divergence is 
cut off at the distances ~ £s where the long-range vortex 
potential <f> (|B19I) vanishes. For b 3> £s we have s c = 0. 
The perturbation approach holds as long as Co <C C and 
C + <C 1, i.e., as long as |6| <C £, 2 d- For s 3> £, 2 d the 
coefficient Co decays faster than exponentially, while 



C+(s)e As -> C_(s)e" 



-As 



2Vr 2 - e 2 P 



such that C is — (6//o)r/-\/r 2 — e 2 and the corrections to 
9 and g vanish as it should be according to (|C1[) . For a 



small distance s — sq (p$ = Sg + b 2 ) we have 

C(so) = C + C + ( So ) + C_( So ), 
0(s o ) - -=L=[r-eC + rc (s ) 



vf 2 " 

+ e[C+(s ) - C_(s )] 



sOo) 



[-e + rC-eC (s ) 



r[c + ( So )-c_( So )]]. 



(CIO) 



(Gil) 



(C12) 



Consider first trajectories that miss the primary vortex 
core, i.e., they go at impact parameters £s <C b <C £, 2 d- 
In this case, the perturbation result Eqs. (|C2IIC4j) can 
be applied along the entire trajectory such that one can 
put so = Sc = 0. The boundary condition for an odd 
function requires 9(0) — 0. Since in this case C_(0) = 0, 
we find from Eq. (|C11|) 

r + eC+(0) =eC-rc* (o) . 

Expressing the coefficients Cq and C+ in terms of the 
energy e 2 (b) of bound states in the induced vortex core, 
Co = -2CC+ = Ce 2 (6)/r, we find 



C[e-e 2 (b)]=T-ee 2 (b)/2T 



(C13) 



where the energy spectrum e 2 (b) of localized excitations 
is given by the Eq. ([5]). According to Eq. (IC13|) e 2 (b) 
is the only spectrum branch in the energy interval |e| <C 
Aoo. The Green function is 



9W 



iTC 



ieC (s) 



Vr 2 



- 
iT 



vT 2 

[c+( 



s e 



As 



vT 2 

C-(s)e 



-Xsl 



(C14) 

For s > £ 2D we have C -> 0, C+e As - C_e~ As -> 0, so 
that the first term is the homogeneous background while 
the rest terms describe the vortex contribution. To ob- 
tain the retarded function for e 2 > T 2 one has to continue 
\/T 2 — e 2 analytically throughout the upper half-plane of 
complex e keeping 3?-\/r 2 — e 2 > 0. 

The normalized LDOS can be found as a sum over 
different trajectories: 



2tv 



N(r,e) 



where s = pcosa', b 



1 

27 



N e (s,b)da' , 
psina', and 



N e ( S) b) = ~(g R (s,b)-g A (s,b)) . 

For |e| < r, a nonzero LDOS comes only from the vor- 
tex contribution of the second and third terms in (|C14[) 
due to the presence of a pole in the coefficient C ac- 
cording to Eq. (|C13|) . The Green functions and LDOS 
reach their long-distance values g = —ie/y/T 2 — e 2 and 
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N = |e|/Ve 2 -r 2 x(e 2 - T 2 ) as p -> oo. For p > £ s 
the trajectories with large impact parameters b > £s 
give the main contribution to the LDOS. In the region 
£s <C P *C £ 2 d we get the angle-resolved density of states 
in the form: 



N e (s,b) 



yr 2 -e 2 (r 2 -e 2 /2) 

fa 

X7T(5[e - e 2 (6)] , 



kl < r (ci5) 

^ (S ' 6) = sign(e)r 2 [ £ - £2 (6)] ' |£| > r (C16) 

Thus, the corresponding LDOS in the energy interval 
|e| < r has the only peaks at e = e 2 (±p): 



1 



db 



7 Vp 2 - b 2 

-p 



vT 2 - e 2 (i - e 2 /2r 2 ) 



X ^(p)- £ 2 ] . (C17) 



For energies above the induced gap, |e| > T, for the same 
distances the LDOS is monotonically increasing with |e| 
to its normal state value: 



7V(p,e) = Ve 2 



r 2 



If! 1 (i- e 2 /2r 2 ) 
2r 2 



00 



(Ci8) 



The LDOS behavior for small distances p < £s de- 
pends crucially on trajectories with small impact param- 
eters b. In this case one has to match Eqs. (|C10|) - (|C12j) 
with the solution obtained in the vortex core region. For 
small s < so we assume the even parts of the Green func- 
tion g and C to be nearly constant, therefore integrating 
Eq. (|B6|) along the trajectory over s from to so we find 
the matching condition for the Green functions: 



hV 2 



2D , 



(*o) = C(*o) 



/■so pso 

/ Ei ds+ig{s Q ) / Eids. (C19) 
Jo Jo 



This matching condition determines the constant C. Its 
poles as a function of energy and the impact parameter 
define the eigenstates of excitations. 

While deriving the effective boundary condition (|C19j) 
f° r b < £s, one needs to separate the exponentially con- 
verging parts Sj°j at s ~ from the long-distance, 
s 3> £s, asymptotics of Eij. For e <C Aoo the long- 
distance expressions, Eq. Q yield Si = 0, E# = Ts/p, 
Ej = Tb/ p. Therefore, we find 



Ei G?S 



E[ oc ds 



(C20) 



/ Eids = I E l j oc ds + T b/pdsK 
Jo Jo Jis 

POO 

« / E^ c d S + rbln( So /£s) . (C21) 
Jo 



The localized self-energy parts Ej°j determine the small- 
distance LDOS and spectrum of excitations. Therefore, 
while Ej°j are dependent on the tunneling model, we 
should consider these models separately. 



1. Coherent Tunneling 

Here we consider the quasiparticle trajectories which 
go through the core of the primary vortex at impact pa- 
rameters b <C £,s assuming coherent tunneling mechanism 
and derive the expressions for the spectrum of localized 
excitations and LDOS. In this case, the self energies are 
equal to the quasiclassical Green functions in the bulk SC 
taken at the same trajectory as in the 2D layer (Fig. 01 : 



-res 



(C22) 



Note that the localized part E 2 oc of the effective or- 
der parameter E 2 has the coordinate dependence E 2 oc = 
iEj (b, s)e ia with zero circulation, unlike its adiabatic 
part gj E 2 (p > £g) = Te^. As we will see below it 
is this different angular dependence of the effective gap 
asymptotics, which leads to the formation of a "shadow" 
of the bulk SC anomalous branch in the excitation spec- 
trum and LDOS in the 2D layer. 

Eqs. (IC7HC9)) yield 



Co(so) = 7 — In- — , 
C + (s )±C.(s ) « — 7— In — ^~ 



ea(&) 

2r 



(C23) 
.(C24) 



We now match the asymptotic solution Eqs. (|C10BC12|) 
obtained for s > s using Eq. ([Cl"9| and Eqs. (jC20l |C2T|) . 
As a result, 



C 



6u[e-e 2 (6)]+2[T- Vr 2 - e 2 



6ur + 2e 



-(r + ^-e 2 ) 



/ Co ds - &D 
Jo 

e 2 (b) 



r 

ee 2 (6) 

2r 

Cods 



Co ds 



(C25) 



r jo 

where Co(s) is the localized part of Cs and 

hvn 

2[e-e (&)] 

Here we put g — iCo and replace the cutoff parameter in 
(O by A = £,2d/£s- For b ^> £s the contributions from 
the primary vortex core proportional to J Q Co ds vanish 
since the trajectory misses the core, and Eq. (|C25|) goes 
over into Eq. ([UTS]) . 

For small b <C £~s the Green function has a pole when 

P(e,b) = [e-e 2 (b)][e-e (b)} 

+ q [r 2 - rVr 2 - e 2 - ee 2 (6)l = (C26) 
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where q = v\\/V2D- This equation coincides with Eq. © 
within the accuracy of our approximation since 62(b) <C 
60(b). The coefficient C takes the form 



C 



[T-ee 2 (b)/2T}[e-e (b)} 
P(e,b) 



g [er- £2 (6)(r + yr^?)/2] 

P(e,b) 



(C27) 



Equation (|C26|) has two real- valued branches of solutions 
61.2(b) in the range |e| < T and one complex branch 61(b) 
in the range T < |e| < A^, for retarded (advanced) Green 
functions. For e <g; T, expanding Eq. (|C26|) in e/r within 
the first order accuracy in 62 (b) we can write 



[6-6 2 (b)}[6~6 a (b)} + p-62(b)} 2 =0 

This equation has two solutions: 

e 1 (6) = (l + (Z /2)- 1 eo(6) 



(C28) 



(C29) 



and 62(b). 

The angle-resolved DOS for small energies |e| <C V and 
p < £s reads 

Nc(s,b) = ^S[6 - 61 (b)] + ^L±^ 5[ 6 ~62 (b)] . 

(C30) 

Here we neglect the terms 662(b) /T 2 and 62(b) /61(b) and 
put eo(b)/ei(6) = l + g/2 according to (|C29|) . In this case 
the LDOS 



N( P ,6) 



Tqx[4(p)-e 2 } T(q + 2) X [el(p) 

^4(p)-t 2 



2V4(P) 



(C31) 



reveals a two-peak structure vs energy at e = £1,2 (p)- For 
|e| ~ r, one can neglect 62(b) and obtain: 



<*(*>)] 



x/r 2 - e 2 



qT6 = 



(C32) 



For |e| > r the dispersion relation is complex valued and 
for retarded functions takes the form: 



e[e - e (b)] + qT V + isign(e)\/ e 2 - T 2 =0. (C33) 

The latter equation describes the resonant states in the 
2D vortex core which decay into the quasiparticle waves 
propagating in the 2D layer above the induced gap. 

Finally, the whole spectrum structure, shown in Fig.JTJ 
has two anomalous branches: (i) one of them 62(6) is 
completely real-valued and follows the CdGM spectrum 
for the superconductor with homogeneous gap T; (ii) an- 
other one is close to the bulk CdGM spectrum, but has 
a discontinuity at e = T, where it becomes essentially 
complex. 



Thus, the LDOS for energies above the induced gap 
|e| > r and small distances p,b < £5 reads 



N(p,6) = 

, <?r 2 



Ve 2 - r 2 



2|e| 



Ve 2 - T 2 - tT 



(C34) 



(e 2 + qT 2 + iqTy/6 2 - T 2 ) 2 - 6 2 6 2 (p) 



and has the only peak at e = 5Rei(p) of the height ~ 
r 2 / e o(/°) f° r P ^ ^sl&D- In the opposite limit of rather 
large distances p > £|/£2D at |e| > T, the spectrum 
reduces to the CdGM spectrum with a finite broadening: 



61(b) = 6 (b) - iTq 



(C35) 



The LDOS has a small difference from its normal state 
value Nq = 1: 



qT 2 

JV(p,e) = l+| ? K 



IT 



^(e + lq r) 2 ~6 2 (p) 



(C36) 



The LDOS in the whole energy range (fC3i~| [C34|l has 
two or even three peaks for such distances. The latter 
case is realized at the distances corresponding to b' < 
b < b* , where the spectrum vs the impact parameter has 
3 anomalous branches. 



2. Incoherent Tunneling 

Assuming small impact parameter values b <C £g, i.e. 
62(b) <C r 2 /A oc , we obtain an expression for the coeffi- 
cient C using the asymptotical solution (|C5[ IC6[) and the 
matching condition (|C19|) : 



C 



6 - 62(b) 



2Vr 2 - e 2 

hv 2D 



Ei ds 



2e 



KV 2 D Jo 



(C37) 



Since |Ei| ~ |Ef c | - T the pole of the coefficient C is 
located at small energies e < T 2 /A <C F. Thus, for e -c T 
the expression for this coefficient takes the form 



C 



£2 (b) 



?2D 



Ei oc ) ds 



r . (C38) 



The localized self energies Ei and Ej oc can be neglected 
for 6 ~ T. They also vanish for |6| 3> £s. In both 
these limits, Eq. (IC37I) transforms into Eq. (|C13|) . The 
integral term in the equation above can be written 
in terms of its real /3(b) = 0i(b) — j3i (b) and imaginary 
7(6) = ji(b) — Ji(b) parts as follows: 



T}j oc ) ds = -/3(b) ± 17(6) . (C39) 



?2D Jo 
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Here upper (lower) sign corresponds to the retarded (ad- 
vanced) Green function. Further we calculate the terms 
of real (3\j and imaginary 71 j parts of the integral (|C39|) . 
which are defined by the following expressions 

00 00 
&»(&) = 7^- / ®X a (s)ds , 7a (6) = / ^ a (s)ds 

t,2D J ?2D J 



and consider the case of the small impact parameter val- 
ues b <C £s : 



2T 2 b f 00 I v\\e~ K 
Jo \ 2Qilp 2 



1 - 



Ve 2 - Wp- 



f/,S 



where p 2 = b 2 + s 2 . In this case the first term in the 
above integral is determined by s ~ 6: 



r& 



00 v\\e- K 



QVLp 2 



ds = Tb 



ds 



o \QO(s 2 + 6 2 ) 

The second one is determined by very small impact pa- 
rameters and reads: 



vw 



IT 

2 



ds 



{s 2 +b 2 )^ 2 ~7 2 ~ 2\be\' 



where fr 2 , = e 2 /fl 2 — b 2 > 0. As a result, we find: 

r 2 /7TUII 

»<») = {qU^ v - ^ 

After simplifying the expression for /3(b) — f3i(b) — /?i(6) 
we obtain (fTTjl . For 6 > £5 the quantity /?(&) decays as 
exp(-26/Cs). 

The expressions for imaginary parts hold for any dis- 
tances p because the delta functions in the integrals select 
only the trajectories that pass at small impact parame- 
ters: 



7i(b) 



2D Jo 



I. 



■,-K 



x (n 2 p 2 -e')) ds 



\Q^/fT 2 p 2 -,- 



A- 



72(6) ^y 



^2D \Qfi V / \fl 2 b 2 -e 2 \ i 

x (n 2 P 2 - e 



wire 



-A" 



/V Qv^V - e 

r 2 ' 



rf.s 



sign(fee) 



«u n|b| 



Here we use the following expressions for the standard 
integrals: 



/' 

Jbn 



ds 



= In- 



A 



/ b0 V^±^ ^/\n 2 b 2 -e 2 \ ' 
where s max — £s, and 



The imaginary terms also decay exponentially for b > 
The expression for 7(6) = 71 (6) — 7/(6) gives (fT2|) . As a 
result, the expression for the coefficient C reads: 

C = r/[e-e 2 (6)-/3(6)+i 7 (6)] (C40) 

and the angle-resolved DOS for e < T takes the form 

T~f(b)e- X ^ 



N £ (s,b) = 



\e-e 2 (b)-f3(b)] 2 +J 2 (b) 



(C41) 



coinciding with Eq. (1101) in the main text. Since param- 
eters /3, 7 ~ r/A and 62(b) /T <C 1 arc small for & -C £2.0 
and |e| > T, the LDOS reaches its bulk value in this 
region: 



N(p,e) 



Ve 2 - r 2 



(C42) 



3. Dirty superconductor and clean 2D Layer 

Here we derive the Green functions and the DOS in 
the 2D layer for the dirty limit of the bulk SC. For small 
impact parameter values b <C £s we get T}° c = and the 
matching condition takes the form: 

00 

& D e(s ) = 2i({s ) j sineds + 2ig(s )bln[s /ts} (C43) 


The coefficient C in this case has the only broadened pole 
at e = £2(6): 



C[e-e 2 {b) + i 1 ]=T , 
where the broadening 



(C44) 



7 



2rvf 



2 _ c 2 r°° 



hv 2D 



sin O ds 



coincides with Eq. (|T5|) in the main text. For |e| < T 
and p < £s the angle-resolved DOS can be written in the 
form 

N e (s,b)= m 2 , ,L M am ■ (C45) 
Vr 2 - e 2 [e - £2(o)J 2 + Yv>) 

Consequently, the LDOS has a peak of the height ~ 
T/j(p) at energy e = e 2 (p). 
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For the energies above the induced gap e > T and small 
impact parameter values e 2 (fe),7(6) -C T the local DOS 
can be replaced by its bulk value: 



N{ Pl e)= W |e| (C46) 

For b 3> £s the imaginary part of energy decays expo- 
nentially, and Eq. (|C44[) transforms into Eq. (|C13|) . 



with £ 2 = S 2 e"^. 

Using a standard parametrization g(p) — T^sin^I + 
t 2 cos ^e~ %T ^ and the expressions for the vortex po- 
tentials one can obtain (|16p from the main text with 
V 2 = p~ 1 d p (pd p ). Integrating Eq. ([To]) , multiplied by 
p, in a small region around the origin (from p — to the 
value < Po < &d) we find the matching condition 
for the adiabatic Green function (|C5[ IC6|) : 



4. Dirty superconductor and dirty 2D layer 

At the end of this section we concentrate our attention 
on the dirty limit both in 2D layer and superconductor. 
In this case the bulk SC (fT5)) and 2D layer Green func- 
tions satisfy the Usadel equations (fT4l HI)]) . Indeed, for 
momentum-orientation-averaged Green functions in 2D 
layer 



S(p) 



.92 

-fie 



-i<p 



fa* 

92 



d 2 p 

{2-Kf 



g(p, r) 



one can derive the equation: 



iD 2D [ 92 (V 2 - p- 2 )h - /2V 2 . 92 ] - 

-2(e+E 1 )/ 2 + 2S 2ff2 =0, (C47) 



D 



2D 



9 T 

op 



sin 2* 



2p 



-dp 



rpo 

2/ pdp[rsin(*-e) + iecos*] =0 . (C48) 
Jo 



Considering the expansion ^(po) = — Kpo with 
K = d 1 i>(po)/dp ~ £^ and assuming ^0 7^ tt/2 one 
obtains cos^o ~ Pa/ivfn m (Po/£,s)) "C 1. This estimate 
confirms the conclusion that the LDOS in the dirty limit 
follows the bulk LDOS pattern scaled with the 2D coher- 
ence length ^2D to within the second order terms in the 
small parameter pof&D- 
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